Here are some examples of the HydroData Package – April 12, 2018.

Examples in this vingette include

  1. USGS streamflow
  2. GHCN
  3. NHD
  4. NLCD
  5. SNOTEL
  6. Weather Underground
  7. WBD
  8. NED
  9. TIGER roads
  10. Airports

NOT shown (but available in package)

  1. Daymet
  2. Waterbodies
  3. Reservoirs
  4. CLD
  5. PRISM
  6. Koppen
  7. NWM
  8. USGS water use
  9. NID

Find USGS Gages near UCSB (in a 20 mile square area…)

 ucsb.usgs = findUSGS(clip_unit = list("ucsb", 20, 20), ids = T, save = T)
## AOI defined as the supplied clip unit.
##            Now loading loading CONUS USGS data...
## All USGS Data loaded: 6 gages in total
## 6 USGS gages found within supplied clip unit
## Returned list includes: USGS NWIS shapefile, and list of station IDs

The save = TRUE parameter saves the shapefile to the user working directory (with respect to the query) for later use in R or Arc:

Explore USGS stations

explore(ucsb.usgs, save = T)

The save = TRUE parameter saves the map to the user working directory for later reference or sharing:

Find GHCN PRCP records near UCSB

ucsb.ppt = findGHCN(clip_unit = list('ucsb', 20, 20), param = "PRCP", ids = T)
## AOI defined as the supplied clip unit. Loading global GHCN data...
## 15 GHCN stations found in supplied clip unit
## Returned list includes: NOAA GHCN shapefile, and list of station IDs

Explore GHCN PRCP stations near UCSB

explore(ucsb.ppt)

Get Flow Records from the USGS stations

(shameless wrapper of the USGS dataRetrival function as indicated in documentation):

ucsb.flow = getUSGS(ucsb.usgs$ids)
head(ucsb.flow)
##   agency_cd  site_no       Date Flow Flow_cd
## 1      USGS 11119745 1983-10-28 0.81       A
## 2      USGS 11119745 1983-10-29 0.60       A
## 3      USGS 11119745 1983-10-30 0.65       A
## 4      USGS 11119745 1983-10-31 0.70       A
## 5      USGS 11119745 1983-11-01 2.20       A
## 6      USGS 11119745 1983-11-02 0.93       A

Get PRCP Records from the first three GHCN stations:

ucsb.ppt = getGHCN(ucsb.ppt$ids[1:3], parameters = "PRCP")
head(ucsb.ppt$PRCP)
##   agency_cd     site_no       Date prcp
## 1      NOAA US1CASB0001 2009-09-01   NA
## 2      NOAA US1CASB0001 2009-09-02   NA
## 3      NOAA US1CASB0001 2009-09-03   NA
## 4      NOAA US1CASB0001 2009-09-04   NA
## 5      NOAA US1CASB0001 2009-09-05    0
## 6      NOAA US1CASB0001 2009-09-06    0

Visualize the Flows

inspect(ucsb.flow, param = "Flow")
## Station 11119745 has 1825 daily FLOW observations.
## Station 11119940 has 1820 daily FLOW observations.
## Station 11120000 has 1824 daily FLOW observations.
## Station 11120500 has 1825 daily FLOW observations.
## Station 11123000 has 1825 daily FLOW observations.
## Station 11123500 has 1825 daily FLOW observations.

Visualize Rainfall

inspect(ucsb.ppt$PRCP, param = "PRCP")
## Station US1CASB0001 has 1825 daily PRCP observations.
## Station US1CASB0004 has 574 daily PRCP observations.
## Station US1CASB0009 has 789 daily PRCP observations.

NHD Example

Here we will access a non-point feature

ucsb.nhd = findNHD(clip_unit = list('ucsb', 20, 20))
## AOI defined as the supplied clip unit. Shapefile determined. Now loading NHD flowline data...
## Trying URL ...
## All NHD flowlines loaded: 422 in total.
## Returned list includes: flowline shapefile

Visualize Flowlines and USGS stations together

explore(list(ucsb.usgs, ucsb.nhd))

Land cover/Raster Example

lc = getNLCD(clip_unit = list("ucsb", 10, 10), year = 2011)
## AOI defined as the supplied clip unit. Shapefile determined. Now loading NLCD data for 2011.
## There is 1 NLCD raster in this scene per year. 1  tiles reuquested...
## Downloading raster 1 of 1
## Finished downloading raster 1 of 1
## Raster number 1 Cropped.

Explore Land cover and USGS stations

explore(list(lc, ucsb.usgs))
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA

Alternative Means of Describing an Area of Interest

By STATE

Find SNOTEL stations in California

ca.snotel = findSnotel(state = 'CA', boundary = T, ids = T)
## AOI defined as the boundary of California. Shapefile determined. Now loading loading Snotel data...
## 34 snotel stations found within boundary of California
## Returned list includes: snotel shapefile, AOI boundary, and list of station IDs
explore(ca.snotel)

Visualize the Snow Water Equivilent for first three stations:

snow.data = getSnotel(ids = ca.snotel$ids[1:3])
## 12263 Snotel values for station 301 downloaded (1/3)
## 25969 Snotel values for station 356 downloaded (2/3)
## 31274 Snotel values for station 1051 downloaded (3/3)
head(snow.data)
##   agency_cd site_no       Date SWE PRCP.ACC TMAX TMIN TAVG PRCP.INC
## 1      NRCS     301 1984-09-13  NA       NA   32   32   32       NA
## 2      NRCS     301 1984-09-14  NA       NA   NA   NA   NA       NA
## 3      NRCS     301 1984-09-15  NA       NA   NA   NA   NA       NA
## 4      NRCS     301 1984-09-16  NA       NA   NA   NA   NA       NA
## 5      NRCS     301 1984-09-17  NA       NA   32   32   32       NA
## 6      NRCS     301 1984-09-18  NA       NA   NA   NA   NA       NA
inspect(snow.data, param = 'SWE')
## Station 301 has 1825 daily SWE observations.
## Station 356 has 1825 daily SWE observations.
## Station 1051 has 1824 daily SWE observations.

By COUNTY

Find Airports in Tuscaloosa County, Alabama

al.ap = findAirports(state = 'AL', county = 'Tuscaloosa', boundary = T)
## AOI defined as the  boundary of Tuscaloosa County, Alabama.
## Shapefile determined.
## Loading Aiport Database
## 1 airports found within  boundary of Tuscaloosa County, Alabama
## Returned list includes: airport shapefile, and AOI boundary
explore(al.ap)

Get WeatherUndergorund Data for airport:

al.data = getWU(airport_code = al.ap$airports$Digit4, year = 2015, month = 1:12, type = 'daily')
## Year 2015 Month 1 downloaded.
## Year 2015 Month 2 downloaded.
## Year 2015 Month 3 downloaded.
## Year 2015 Month 4 downloaded.
## Year 2015 Month 5 downloaded.
## Year 2015 Month 6 downloaded.
## Year 2015 Month 7 downloaded.
## Year 2015 Month 8 downloaded.
## Year 2015 Month 9 downloaded.
## Year 2015 Month 10 downloaded.
## Year 2015 Month 11 downloaded.
## Year 2015 Month 12 downloaded.
head(al.data)
##         Date Year Month Day T_max T_avg T_min DP_max DP_avg DP_min H_max
## 1 2015-01-01 2015     1   1    54    43    31     42     35     29    92
## 2 2015-01-02 2015     1   2    55    50    44     53     47     42   100
## 3 2015-01-03 2015     1   3    68    62    55     66     62     54    93
## 4 2015-01-04 2015     1   4    65    54    42     63     46     32   100
## 5 2015-01-05 2015     1   5    49    40    30     31     25     22    85
## 6 2015-01-06 2015     1   6    56    42    27     33     28     25    92
##   H_avg H_min SLP_max SLP_avg SLP_min V_max V_avg V_min Wind_max Wind_avg
## 1    69    46   30.42   30.34   30.27    10     9     5        8        3
## 2    93    86   30.34   30.23   30.14    10     3     1       12        5
## 3    90    87   30.12   30.00   29.92    10     5     1       17        8
## 4    77    54   30.51   30.18   29.97    10    10     4       20       11
## 5    60    34   30.68   30.56   30.42    10    10     8       17        8
## 6    64    35   30.47   30.40   30.32    10    10     7       14        4
##   Wind_min PPT_tot                  Events          agency_cd site_no
## 1       10    0.30                    Rain WeatherUnderground    KTCL
## 2       14    1.61 Rain\n\t,\nThunderstorm WeatherUnderground    KTCL
## 3       23    2.47 Rain\n\t,\nThunderstorm WeatherUnderground    KTCL
## 4       28    0.37                    Rain WeatherUnderground    KTCL
## 5       23    0.00                         WeatherUnderground    KTCL
## 6       19    0.00                         WeatherUnderground    KTCL

Visualize average daily wind velocity

inspect(al.data, param = "Wind_avg" )
## Station KTCL has 365 daily WIND_AVG observations.

Change timesetep to monthly

inspect(al.data, param = "Wind_avg", timestep = 'monthly' )
## Station KTCL has 12 monthly WIND_AVG observations.

By a user-supplied shapefile (here all Los Angeles communities subset to Santa Monica)

LA = rgdal::readOGR("/Users/mikejohnson/Documents/GitHub/Communities1/Communities.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/mikejohnson/Documents/GitHub/Communities1/Communities.shp", layer: "Communities"
## with 353 features
## It has 13 fields
## Integer64 fields read as strings:  X_CENTER Y_CENTER
SM = LA[LA$LABEL_CITY == 'Santa Monica',] %>% spTransform(HydroDataProj)

leaflet() %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons(data = SM,fillColor = "transparent", color = "black", weight = 4)

Find all HUC 10 units intersecting Santa Monica

watersheds = findWS(clip_unit = SM, level = 10)
## Determining intersecting HUC8 units...
## There are 1 HUC8 units in this AOI: 18070104
## 
explore(watersheds)

Find Elevation data for the Ballona Creek Watershed (HUC10)

sm.elev = getNED(clip_unit = watersheds$huc_10[1,], res = 1)
## AOI defined as the supplied clip unit. Shapefile determined. Now loading 1 arc second NED data...
## There are 2 NED rasters in this scene.
## Downloading raster 1 of 2
## Finished downloading raster 1 of 2
## Downloading raster 2 of 2
## Finished downloading raster 2 of 2
## Raster number 1 Cropped.
## Raster number 2 Cropped.
## Mosaicing raster...
## Returned object contains elevation raster

Plot and ‘analyze’ elevation data for the Ballona Creek Watershed (HUC10)

explore(sm.elev)
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
hist(sm.elev$elev, xlab = 'Elevation (m)', main = paste("Elevation Distribution", watersheds$huc_10[1,]$NAME, "Watershed"), col = 'black')

Find the roads of Santa Monica

sm.roads = findRoads(clip_unit = list("Santa Monica", 2, 2))
## AOI defined as the supplied clip unit. Loading TIGER roads database...
## There is 1 TIGER file in this scene.
## Trying URL ...
## All TIGER 1 loaded: 137,781 in total.
## Returned list includes: TIGER roads shapefile
explore(sm.roads)